*WIC co-op eWIC TWFE diff in diff analysis
*last modified: 25 August 2023
*last modified by: Charlotte Ambrozek
*this do file implements a standard twfe regression to estimate effects of WIC EBT on all outcomes

clear all
set more off
set rmsg on

set maxvar 120000
set emptycells drop

local data_dir ./data/cleaned
local raw_dir ./data/raw
local out_dir ./analysis/output
local graph_dir ./analysis/output/graphs
local tab_dir ./analysis/output/tables
local log_dir ./documentation/logs
local date: display %tdYY-NN-DD date(c(current_date), "DMY")
di "`date'"
capture log close

log using `log_dir'/ewic_twfe_all`date', replace

*import cleaned eWIC implementation info into own frame
*eWIC implementation variable is at county fips/fiscal year level 
*partition data into state/treatment year groups

*import cleaned ewic implementation info into own frame
mkf ewic_fy
frame ewic_fy{
	use ./data/cleaned/ewic_rollout.dta, clear
	rename year fiscalyear
	drop pre* post* 
	*fill in gaps in state fips using modes
	generate str_fips = string(ct_fips)
	replace str_fips = "0" + str_fips if strlen(str_fips) == 4
	generate st_fips = substr(str_fips, 1, 2)
	destring st_fips, replace 
	bys st_fips: egen state_mode = mode(state)
	replace state = state_mode if state != state_mode & !missing(state_mode)
	drop str_fips state_mode

	*drop Vermont, Mississippi, Nevada,  + missing st_fips
	drop if st_fips == 50 | st_fips == 32 | st_fips == 28 | missing(st_fips)
	duplicates report ct_fips fiscalyear 
	assert r(N) == r(unique_value)
	compress
	frame put ct_fips fiscalyear ev_year state st_fips, into(ct_sq_frame)
}

*import TIP store data into own frame
mkf tip_fy
frame tip_fy{
	local data_dir ./data/cleaned
	use `data_dir'/tip_auth_sq, clear
	*drop "Direct Distribution Center" and "Home Food Delivery Contractor" types
	gen f = (vendor_type1 == 3 & (vendor_type2 == 4 | missing(vendor_type2))) | (vendor_type1 == 4 & (vendor_type2 == 3 | missing(vendor_type2)))
	bys tip_id: egen flag = mode(f)
	bys tip_id: replace vendor_type1 = vendor_type1[_n+1] if flag == 0 & f == 1
	drop if flag == 1
	drop flag f tip_year_id
	*fill in gaps in state fips using modes
	generate str_fips = string(ct_fips)
	replace str_fips = "0" + str_fips if strlen(str_fips) == 4
	generate st_fips = substr(str_fips, 1, 2)
	destring st_fips, replace 
	
	*drop Vermont and mississippi because of different regiemes prior to eWIC
	*drop NV because ewic turn on then off then on again and we don't have reliable timings
	*drop missing st_fips
	drop if st_fips == 50 | st_fips == 32 | st_fips == 28 | missing(st_fips)

	*flag stores that specialize in WIC
	gen a50 = vendor_type1 == 1 | vendor_type1 == 7 | vendor_type2 == 1 | vendor_type2 == 7

	*link ewic implementation info
	frlink m:1 ct_fips fiscalyear, frame(ewic_fy) generate(ewic_link)
	frget ev_year , from(ewic_link)
	* assume treated at earliest exposure
	bys tip_id: egen ey_min = min(ev_year)
	* assume treated at earliest exposure
	replace ev_year = ey_min if ev_year != ey_min & !missing(ey_min)
	bys tip_id (ev_year): assert ev_year[1] == ev_year[_N]

	gen ebt = fiscalyear >= ev_year

	*assert no dups at future level of xtset 
	duplicates report tip_id fiscalyear 
	assert r(N) == r(unique_value)
	xtset tip_id fiscalyear
	compress

	*this zip has multiple implementation dates and can't be reconciled
	drop if zip == 42223
	
	*put fields necessary to create a frame of unique zip/year obs with ewic implementation; later will merge this to tip_fy again and then collapse to get aggregated zip level outcomes
	compress
	frame put zip fiscalyear ev_year state st_fips, into(zip_sq_frame)
	*put all obs into frame to use for zip collapsing
	frame put _all, into(zip_sq)
}

*import two separate cleaned redemptions datasets (second is balanced to only include ZIPs that have non-missing redemptions in all years)
mkf red_fy
frame red_fy{
	use ./data/cleaned/tip_wic_redemptions_05_18, clear
	rename (wic_redemptions) (wicred)
}

mkf red_bal_fy
*load balanced redemptions as an attempted fix
frame red_bal_fy{
	use ./data/cleaned/tip_wic_bal_redemptions_05_18, clear
	rename (wic_redemptions) (wicred)
}

frame zip_sq_frame{
	duplicates drop 
	
	bys zip: assert state[1] == state[_N]
	bys zip: assert st_fips[1] == st_fips[_N]
	*resolve discrepancies in values of ewic implementation variables as follows:
	*ev_year is equal to the min of ev_year across the zip/fiscal year (use the earliest possible event year)
	bys zip: egen ey_min = min(ev_year)
	* assume treated at earliest exposure
	replace ev_year = ey_min if ev_year != ey_min & !missing(ey_min)
	bys zip: assert ev_year[1] == ev_year[_N]

	duplicates drop

	duplicates report zip fiscalyear 
	assert r(N) == r(unique_value)	
	compress
}

*collapse to zip/fiscalyear level and then link ewic policy variables back in
frame zip_sq{
	*reshape by chain status to allow looking at chain status as heterogeneity later; may also want to do this for a50 or vendor type at some point
	collapse (sum) auth (first) state st_fips, by(zip fiscalyear chain)
	gen cstring = ""
	replace cstring = "chain" if chain == 1
	replace cstring = "indep" if chain == 0
	drop chain
	reshape wide auth, i(zip fiscalyear) j(cstring) string
	assert !missing(zip) & !missing(fiscalyear)
	egen auth = rowtotal(auth*)
	recode authindep authchain (missing = 0)
	*start linking in other policy variables and outcomes (redemptions)
	frlink 1:1 zip fiscalyear, frame(zip_sq_frame) generate(ewic_link)
	frlink 1:1 zip fiscalyear, frame(red_fy) generate(red_link)
	frlink 1:1 zip fiscalyear, frame(red_bal_fy) generate(red_bal_link)
	frget ev_year , from(ewic_link)
	frget wicred, from(red_link)
	frget wicbalred = wicred, from(red_bal_link)
	compress
	assert !missing(ev_year) & !missing(fiscalyear) & !missing(zip)
	xtset zip fiscalyear
    gen ebt = fiscalyear >= ev_year
	compress
	*log 
	foreach var of varlist auth authindep authchain {
		format `var' %11.6f
		local r = abs(rnormal(0, 0.000001))
		tempvar zero
		gen `zero' = `var' == 0
		replace `var' = `r' if `zero' == 1
		gen double l`var' = log(`var')
		*check that log conversion using errors worked
		assert !missing(l`var') if !missing(`var')
		replace `var' = 0 if `zero' == 1
 	}
	foreach var of varlist wicred wicbalred{
		gen double l`var' = log(`var')
		bys zip: egen std`var' = std(`var')
		gen double `var'hk = `var'/100000
	}
	*standardize vars to see if this helps with left censoring
	foreach var of varlist auth authindep authchain{
		*use built in egen fn for this (mean zero, sd 1)
		bys zip: egen std`var' = std(`var'), mean(0) sd(1)
	}
	compress
}

*import cleaned SNAP authorized stores from STARS data
mkf snap_stores
frame snap_stores{
	use `data_dir'/stars_stores, clear
	cap drop __00*
	*this zip has multiple implementation dates and can't be reconciled
	destring zip, replace
	drop if zip == 42223
	compress

	codebook st state store_type

	drop if st == "VT" | st == "MS" | st == "NV"

	*assert no dups at future level of xtset 
	duplicates report stars_id fiscalyear 
	assert r(N) == r(unique_value)
	*this zip has multiple states and can't be reconciled
	drop if zip == 42223 
	*start linking in other policy variables and outcomes (redemptions)
	frlink m:1 zip fiscalyear, frame(zip_sq_frame) generate(ewic_link)
	frget ev_year st_fips, from(ewic_link)
	
	compress
	drop if missing(ev_year)
	*impose balance on set of stars stores (need to know their auth history for full sample)
	bys stars_id: egen yrs = count(fiscalyear)
	assert yrs == 14
	drop yrs
	assert !missing(ev_year) & !missing(fiscalyear) & !missing(st_fips)
	xtset stars_id fiscalyear
	compress
}

*import two separate cleaned redemptions datasets (second is balanced to only include FIPS that have non-missing redemptions in all years)
mkf snap_red_fy
frame snap_red_fy{
	use ./data/cleaned/ct_mon_redemptions, clear
	rename fips5 ct_fips
	collapse (sum) snapred = redemptions, by(ct_fips fiscalyear)
}	

mkf snap_red_bal_fy
frame snap_red_bal_fy{
	use ./data/cleaned/ct_mon_balanced_redemptions, clear
	rename fips5 ct_fips
	collapse (sum) snapbalred = redemptions, by(ct_fips fiscalyear)
}

frame ct_sq_frame{
	duplicates drop 
	
	bys ct_fips (state): assert state[1] == state[_N]
	bys ct_fips (st_fips): assert st_fips[1] == st_fips[_N]
	bys ct_fips: assert st_fips == floor(ct_fips/1000)
	*resolve discrepancies in values of ewic implementation variables as follows:
	*ev_year is equal to the min of ev_year across the ct_fips/fiscal year (use the earliest possible event year)
	bys ct_fips: egen ey_min = min(ev_year)
	* assume treated at earliest exposure
	replace ev_year = ey_min if ev_year != ey_min & !missing(ey_min)
	drop ey_min
	bys ct_fips (ev_year): assert ev_year[1] == ev_year[_N]

	*drop Vermont, Mississippi, Nevada, + missing st_fips
	drop if st_fips == 50 | st_fips == 32 | st_fips == 28 | missing(st_fips)

	duplicates drop
	assert !missing(ev_year) & !missing(fiscalyear)

    gen ebt = fiscalyear >= ev_year

	duplicates report ct_fips fiscalyear 
	assert r(N) == r(unique_value)	

	*start linking policy variables and outcomes (redemptions)
	frlink 1:1 ct_fips fiscalyear, frame(snap_red_fy) generate(red_link)
	frget snapred, from(red_link)
	frlink 1:1 ct_fips fiscalyear, frame(snap_red_bal_fy) generate(red_bal_link)
	frget snapbalred, from(red_bal_link)
	gen double lsnapred = log(snapred)
	bys ct_fips: egen stdsnapred = std(snapred), mean(0) sd(1)
	gen double snapredmil = snapred/1000000
	gen double lsnapbalred = log(snapbalred)
	bys ct_fips: egen stdsnapbalred = std(snapbalred), mean(0) sd(1)
	gen double snapbalredmil = snapbalred/1000000
	xtset ct_fips fiscalyear
	drop if missing(ev_year)
	compress
}

cwf snap_stores
keep stars_id fiscalyear ev_year st_fips auth chain store_type
gen ebt = fiscalyear >= ev_year

eststo snapall: areg auth ebt i.fiscalyear, absorb(stars_id) vce(cluster st_fips)

eststo snapchain: areg auth ebt i.fiscalyear if chain == 1, absorb(stars_id) vce(cluster st_fips)

eststo snapindep: areg auth ebt i.fiscalyear if chain == 0, absorb(stars_id) vce(cluster st_fips)

cwf zip_sq
***fixed effects will be zip code and year, regressing wicredhk on ebt indicator
eststo wicred: areg wicredhk ebt i.fiscalyear, absorb(zip) vce(cluster zip)

cwf tip_fy
frame drop ewic_fy

cwf tip_fy
cd `tab_dir'
eststo tipall: areg auth ebt i.fiscalyear, absorb(tip_id) vce(cluster st_fips)

eststo tipchain: areg auth ebt i.fiscalyear if chain == 1, absorb(tip_id) vce(cluster st_fips)

eststo tipindep: areg auth ebt i.fiscalyear if chain == 0, absorb(tip_id) vce(cluster st_fips)

cwf ct_sq_frame

***fixed effects will be county fips and year, regressing snapredmil on ebt indicator
eststo snapred: areg snapredmil ebt i.fiscalyear, absorb(ct_fips) vce(cluster ct_fips)

esttab tipall tipchain tipindep snapall snapchain snapindep using ewic_twfe_auth_tab.tex, replace tex mtitles("All" "Chain" "Indep." "All" "Chain" "Indep.") numbers keep(ebt) indicate("FY FE = *.fiscalyear") b(3) ci(3) coeflabels(ebt "WIC EBT") mgroups("WIC authorization" "SNAP authorization", pattern(1 0 0 1 0 0)) 

esttab wicred snapred using ewic_twfe_red_tab.tex, replace tex mtitles("WIC" "SNAP") numbers keep(ebt) indicate("FY FE = *.fiscalyear") b(3) ci(3) coeflabels(ebt "WIC EBT") mgroups("Redemptions", pattern(1 0)) 

log close
